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^ Abstract 

illustrate through exphcit numerical calculations how the Born-rule probability densities of non-relativistic quan- 

i-^um mechanics emerge naturally from the particle dynamics of de Broglie-Bohm pilot-wave theory. The time evo- 
Q-i 

-4-4ution of a particle distribution initially not equal to the absolute square of the wave function is calculated for a 

particle in a two-dimensional infinite potential square well. Under the de Broglie-Bohm ontology, the box contains 

i__an objectively-existing 'pilot wave' which guides the electron trajectory, and this is represented mathematically by a 

CNBchrodinger wave function composed of a finite out-of-phase superposition of M energy eigenstates (with M ranging 
> 

O^rom 4 to 64). The electron density distributions are found to evolve naturally into the Born- rule ones and stay there; 

^n analogy with the classical case this represents a decay to 'quantum equilibrium'. The proximity to equilibrium 

„.e. , ... . ... . _ .... 

^ towards zero over the course of time. The timescale r for this relaxation is calculated for various values of M and the 
^oarse-graining length e. Its dependence on M is found to disagree with an earlier theoretical prediction. A power 

^aw, r oc M~^, is found to be fairly robust for all coarse-graining lengths and, although a weak dependence of r on £ 
is observed, it does not appear to follow any straightforward scaling. A theoretical analysis is presented to explain 
these results. This improvement in our understanding of timescales for relaxation to quantum equilibrium is likely 
to be of use in the development of models of relaxation in the early universe, with a view to constraining possible 
violations of the Born rule in inflationary cosmology. 



I. INTRODUCTION 



The Born rule is the fundamental connection between the mathematical formalism of quantum theory and the 
results of experiments. It states that if an observable corresponding to a Hermitian operator A is measured 
in a system with pure quantum state the probability of an eigenvalue Aj will equal ^\E'|-Pj|\l/^, where 
Pi is the projection onto the eigenspace of A corresponding to A^. For the case of measurement of the 
position X of - say - an electron in a box, the probability density at time t for finding the electron at x 
is p{'x,t) = |\l/(x, The Born rule is normally presented as a postulate, though attempts to derive it 
from more fundamental principles have a long history. There has, for example, been much recent work on 
deriving the Born rule within the framework of the many-worlds interpretation of quantum mechanics, but 
such derivations remain controversial [Ij. According to a recently-published (2008) encyclopedia of quantum 
mechanics "t/ie conclusion seems to be that no generally accepted derivation of the Born rule has been given 
to date, but this does not imply that such a derivation is impossible in principle" [2]. 

Born's 1926 paper [3], and Heisenberg's introduction of the uncertainty relations the following year [4j, were 
instrumental in popularizing the idea that Nature at the quantum level is fundamentally probabilistic. The 
idea that \1/ provides a complete description of a single electron strongly suggests that the probabihstic 
interpretation of expresses an irreducible uncertainty in electron behaviour that is intrinsic in Nature. 
It is somewhat ironic therefore - and unknown to most physicists - that Born's rule emerges quite naturally 
out of the dynamics of a deterministic process that was first outlined by de Broglie in 1927 [5]. The 
process in question can be described by a theory commonly referred to as the de Broglie-Bohm 'pilot-wave' 
formulation of quantum mechanics While the theory has attracted little serious interest ever since 

it was introduced [21 E], there has been a considerable resurgence of activity in this area over the last 
fifteen years or so [15]. One of the reasons for this is that, although the theory is completely consistent with 
the full range of predictive-observational data in quantum mechanics, it also permits violations of the Born 
rule and, at least in principle, this leads to the possibihty of new physics and of experimentally testable 
consequences [T6HT9] . 

The reason that de Broglie-Bohm theory can get away with such an apparently absurd contradiction of one 
of the basic postulates of quantum theory is that it assumes orthodox quantum mechanics is incomplete, 
as Einstein always insisted. It supposes that electrons, for example, are real 'particles' with continuous 
trajectories and that the Schrodinger wave function represents an objectively-existing 'pilot wave' which 
turns out to influence the motion of the particles. Since the particle density p and the square of the pilot 
wave are logically distinct entities, they can no longer be postulated to be equal to each other. Rather, 



their identity should be seen as dynamically generated in the same sense that one usually regards thermal 
equilibrium as arising from a process of relaxation based on some underlying dynamics (though with a 
dynamics on configuration space rather than phase space). Since pilot-wave theory features a different set of 
basic axioms and conceptual structures, with event-by-event causality and the prospect of making predictions 
different from orthodox QM, it is better to think of it as a different theory, rather than a mere 'interpretation' 
of quantum mechanics. 

In the pilot-wave formulation then, quantum mechanics emerges as the statistical mechanics of the underlying 
deterministic theory. If the particle distribution obeys the Born rule p = |\E'| , the system is said to be in 
'quantum equilibrium'. One finds in general that: 

1. Non-equilibrium systems naturally tend to become Born-distributed over the course of time, on a coarse- 
grained level, provided the initial conditions have no fine-grained microstructure |16l I20H22] . (The latter 
restriction is similar to that required in classical statistical mechanics. An assumption about initial conditions 
is of course required in any time-reversal invariant theory in order to demonstrate relaxation [211 [22].) 



2. Once in quantum equilibrium a system will remain in equilibrium thereafter, as was originally noted by 
de Broglie [6] (this property is sometimes referred to as 'equivariance'). 

These two observations - along with a description of how these statements about the objective makeup of 
the system might be translated into statements about measurement - can be said to 'explain' or derive the 
Born rule. Given the common general viewpoint referred to in the first paragraph, many physicists might 
consider this surprising. 

In this work, we present a numerical analysis of the timescale for the relaxation of non-equilibrium distribu- 
tions of particles to Born-rule quantum equilibrium using pilot-wave dynamics; the approach to equilibrium 



is monitored by computing the coarse-grained subquantum i/-function (see section IB and Ref. [T6|) . The 
results we obtain are for a particle in a 2D infinite potential square well where the wave function is a finite 
superposition of M eigenfunctions (where, depending on the choice of initial state, M ranges from 4 to 64). 
The initial particle distribution is deliberately chosen to be 'out of equilibrium' by giving it the same form as 
the absolute square of the ground-state wave function, that is, p = 4/7r^ sin^ x sin^ y. This system - with fixed 
M - has been studied before in this context by Valentini and Westman [22j and by Colin and Struyve [23] , 
but here we go further. Our recent development of a new and much faster computer code |24j allows us to 
study systems with many more modes. The timescale r for relaxation is studied as a function of the number 
of modes M (and, in consequence, of the number of nodal points in the wave function) and as a function of 
the coarse-graining length e. The dependence of the relaxation timescale on these two quantities is compared 



to theoretical predictions. 

It is intended that calculations such as these will provide a next step towards a detailed understanding of 
relaxation to quantum equilibrium in the early universe, with a view to constraining possible non-equilibrium 
effects in cosmology. In Ref. [19] it was shown, in the context of inflationary cosmology, that corrections to 
the Born rule in the early universe would in general have potentially observable consequences for the cosmic 
microwave background (CMB). This is because, according to inflationary theory, the primordial perturbations 
that are currently imprinted on the CMB were generated at early times by quantum vacuum fluctuations 
whose spectrum is conventionally determined by the Born rule. To make detailed predictions for possible 
anomalies in the CMB, however, requires a precise understanding of how fast relaxation would occur in, for 
example, a pre-inflationary era (as discussed in section IV-A of Ref. [19]) . It may be hoped that numerical 
studies, such as those reported in this paper, will reveal how the relaxation timescale depends on general 
features of the quantum state such as the number M of modes in a superposition. The results could then be 
applied in future work to specific cosmological models. 



A. Pilot-wave dynamics 



The basic ideas of de Broglie-Bohm pilot-wave theory may be simply understood in a non-relativistic con- 
text [25j. It is a non-local hidden- variables theory, that is, the theory contains some variables that distinguish 
the individual members of an ensemble that in orthodox QM would be considered identical since they all have 
the same wave function. These variables are supposed to be ultimately responsible for the apparently random 
nature of - for example - position measurements on the system. If, as required by some interpretations, one 
were to suppose both that a complete description of the system is afforded by \I' and that \1/ has an objective, 
physical existence, one might conclude from the results of measurements that Nature is intrinsically prob- 
abilistic or random. In pilot-wave theory, by contrast, one supplements the wave function description with 
'hidden variables' by postulating the existence of particles with definite positions, in addition to the wave. 
These particles then follow deterministic trajectories (the nature of which can be deduced) and the observed 
randomness is then understood to be a consequence merely of our ignorance of the initial conditions, that is, 
the starting positions of the particles. 

How does an individual quantum system evolve in time? The pilot wave evolves at all times according to 



the usual time-dependent Schrodinger equation 



i=l 

As normally understood the evolving quantum system behaves like a 'probability fluid' of density = 
with an associated time-dependent quantum probability current, deflned in the usual manner as j = 
^Im(^*V^). In pilot-wave theory, the particles have a continuous objective existence, with trajectories that 
follow the streamlines of the current. Thus their velocity is given by the current divided by the density, that 
is, by 

h 

V = — ImVln\l/. 

m 

Using the complex polar form of the wave function = |\E'| exp[z5'//i], we can recover the (locally deflned) 
phase ^(xi, . . . ,XAr,t) of the wave by the expression S = hlmln"^. The de Broglie guidance equation for the 
trajectories Xj(t) may then be written as 

dt ITLi 

If, for an ensemble of particles with the same wave function, the initial positions have a Born-rule distribution, 
then (by construction) the law of motion of Eqn. [l] implies that the particle positions will have a Born-rule 
distribution at all times. 

If desired, one may take the flrst time derivative to write the equation of motion in second-order form, 

m^x, = -Vi(y + Q), (2) 

where the quantum potential Q = — Xlj2m"n*F' ^^^^ approach, the system acts as if there were a 
'quantum force' —ViQ acting on the particles in addition to the classical force —ViV. This second-order 
approach with a law of motion given by Eqn. [2] was proposed by Bohm in 1952. It may be referred to 
as 'Bohm's dynamics' in order to distinguish it from 'de Broglie's dynamics' based on Eqn. [l] (which was 
proposed by de Broglie in 1927). For de Broglie, p = VS is the law of motion; for Bohm - at least the Bohm 
who wrote the 1952 papers - it is an initial condition which can be dispensed with (clearly if we integrate 
the second-order formula we only recover de Broglie's equation up to some constant and this must be flxed 
for each trajectory by some boundary condition, such as that implied by de Broglie's equation for some time 
to). Thus, in principle, Bohm's dynamics encompasses what one might call 'extended nonequilibrium' where 



p 7^ VS in addition to p 7^ Recent work [26] suggests that this 'extended nonequihbrium' is unstable 
and does not relax in general; if this is correct then it may be argued that Bohm's second-order dynamics is 
untenable as a fundamental theory as there would be no reason to expect equilibrium in the universe today, 
and that de Broglie's dynamics is in fact the fundamental formulation of pilot-wave theory. 

Some additional relevant observations: 

(1) The form of the guidance equation may be altered, while retaining consistency with the Born-rule distri- 
bution. This can be achieved by adding a divergence-free term (divided by to the right-hand side. Such 
alternative velocity fields will not be discussed further here but have been studied by, for example, Colin 
and Struyve [23J and Timko and Vrscay [2Z]. Note that such alternatives yield an equivalent physics only 
in the equilibrium state; away from equilibrium, 'subquantum' measurements would allow one to track the 
trajectories and so distinguish the true velocity field [TTj. 

(2) Given the wave function for a system, the particle trajectories from any starting point may be calculated 
using only the initial position of the particle, rather than the position and the momentum. This is because the 
guidance equation alone gives the particle velocity and consequently the momentum for any initial position. 

(3) Particle trajectories tend to be quite erratic, even with simple wave functions that are superpositions of 
just a few energy eigenfunctions. Fig. [T] illustrates the divergence of neighbouring particle trajectories by 
showing the paths of two particles with almost identical initial positions, propagating according to pilot-wave 
dynamics. 

How do numerical simulations demonstrating the Born rule for the actual particle positions translate into 
statements about 'measurement'? Ideal measurements of position in pilot-wave theory are usually correct 
measurements (they reveal the pre-existing position of the particle - see Ref. [9], p. 351) and so the Born 
rule in position space follows immediately if the particles really are distributed that way. For other kinds of 
measurements, a clear derivation of the Born rule may be found in Section 8.3.5 of Holland's textbook [9] 
(noting that Holland assumes the distribution of actual particle positions as a postulate - see Ref. HJ 
p. 67). The key point is that, in a theory of particles, experimental observations may be reduced to particle 
positions (dots on screens, apparatus pointer positions, etc.) - where laboratory apparatus is treated as just 
another system made of particles. As long as the Born rule holds for the joint distribution of positions of all 
the particles involved (including the particles making up the equipment), then the marginal probability dis- 
tribution for, say, pointer positions (obtained by integrating out the other degrees of freedom) will necessarily 
match the predictions of quantum mechanics. In such a case, the distribution of macroscopically-recorded 
outcomes will be the same as in quantum theory. 
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FIG. 1. Comparison of two trajectories with almost identical initial positions, shown individually in 1(a) and|l(b) for 



clarity, then superimposed in l(c)[ in a system where the wave function is a superposition of 16 energy eigenfunctions. 



The particles have been propagated for time An in both cases ~ note the rapid divergence. |l(a) also shows a good 



example of particle motion near to a node. The trajectory is seen to spiral around a moving nodal point before 
departing from the vicinity of the node (similar behaviour is reported in Ref. I22l) . This behaviour seems to be a 
major driving force behind relaxation. 



B. Quantum equilibrium 

To demonstrate equivalence to quantum mechanics, it is usually simply assumed that the actual distribution 
of particle positions is already supplied to us obeying the Born rule p = |\[^| . In the approach taken here, 
where we try to demonstrate why this is so, the Born-rule distribution is considered to be a special case and 



the particles are said to be in quantum equilibrium when in this state. The dynamics described in section TX 
can just as well be used to describe the evolution of non-equilibrium systems, whereas standard formulations 
cannot. In general in such studies, the probability density is found to approach the Born rule distribution 
over time; it is said to relax to equilibrium [161 ESj. This relaxation is a consequence of the deterministic 
motion of the particles and is not an intrinsically stochastic process (further insight into relaxation has been 
obtained by Bennett using techniques from Lagrangian fluid dynamics [28]). 

Fig. [2] shows the results of a numerical simulation of this relaxation process. It can be clearly seen that 
the particle distribution p rapidly comes to resemble the (periodically repeating) time-dependent The 
example chosen - a superposition of sixteen modes, for a particle moving in two spatial dimensions - is identical 
to that studied by Valentin! and Westman |22]. The results obtained match theirs, thereby providing an 



important confirmation of tlie previous results, witli an independently written and implemented numerical 
code [22] ■ 
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(a)t=0 (b)t=27r (c)t=47r 

FIG. 2. Computational model of the relaxation of an initially non-equilibrium distribution, p = sin^ x sin^ y, 
evolving according to pilot-wave dynamics. The wave function is a superposition of the first sixteen eigenstates for 
a particle in a 2D infinite potential square well. The simulation was run for one period of the wavefunction, or 47r 
in these units. Even after such a short period, significant relaxation towards equilibrium can be observed. (These 
results provide an independent confirmation of those first obtained by Valentini and Westman |22j.) 



If we are to have any chance of observing new physics associated with quantum nonequilibrium states [TBl - 
|2T| [30] , it is important to have a handle on the timescale of this relaxation. To quantify the proximity of a 
distribution to equilibrium, we may use an analogue of the classical //-function [161 122]. This 'subquantum 
//-function' is defined as 




where the (weighted) volume element dS = I^E'I^ dx. 

This quantity will be zero if and only if p = |\E'|^ everywhere, and will be positive otherwise [31], making it 
a useful measure of proximity to equilibrium. Clearly, H is simply the negative of the relative entropy of p 
with respect to |\E'| . 

A feature of the quantities in this definition is that both the volume element, dE, and the ratio / = are 
preserved along trajectories. To show this for /, consider the two continuity equations: 

^ + V-(ip) = (4) 
which follows from the assumption that the actual trajectories follow the velocity field given by Eq. [T], and 

V-(i|^|^)=0 (5) 



dt 

which follows from the Schrodinger equation. 

These two equations can be used to show that the ratio / = obeys: 

Thus / will be preserved along trajectories. Thus, if the system is initially in quantum equilibrium, with 
/ = 1 everywhere, it will never depart from that state. This can, of course, be seen directly from the fact that 
p and |\l'p obey identical continuity equations: if p and |\E'p are initially equal, they will necessarily remain 
equal at all times, since their time evolutions are determined by the same partial differential equation. 



For general (non-equilibrium) initial conditions, the exact value of the i7-function remains unchanged as the 

I 1 2 I 1 2 "i i~2^ 

system evolves. However, if a coarse-graining is applied to p and |\E'| , that is, we replace p — )■ p, |^| — J- I^E'I 
where overbars indicate averaging over small coarse-graining cells, then the coarse-grained i/-function 



i7=/pln(^jdx (7) 

can be shown to be non-increasing, on the assumption that the initial state contains no fine-grained mi- 
crostructure (as in the analogous classical coarse-graining i7-theorem) |26j. Furthermore, H will in fact 
decrease, if the initial velocity field varies with position across the coarse-graining cells [20l|25. The decrease 
of H represents a relaxation of the system towards equilibrium, and formalizes an analogue of the intuitive 
idea of Gibbs: an initial non-equilibrium distribution will tend to develop fine-grained microstructure and 
become closer to equilibrium on a coarse-grained level. Heuristically speaking, this may be thought of in 



terms of two 'fluids', with densities p and that are 'stirred' by the same velocity fleld, and thereby tend 
to become indistinguishable when coarse-graining is applied. 

The effects of coarse-graining on the particle density at some randomly-selected time may be seen in Fig. |3j 








(a)fine- grained 



(b)coarse-grained, e = 32 



(c)smoothed 






(d)e = 32 (e)single coarse-graining cell (f)smoothing applied 

FIG. 3. Figures illustrating the effects of coarse-graining for a 16-mode system at t = ^. 



3 a 



shows a snapshot of 



the fine-grained density p on a 1024 x 1024 lattice; |3(b) shows the coarse-grained density derived from averaging 
the fine-grained density over square cells containing 32 x 32 lattice points (in which case we say the system has a 



'coarse-graining length e = 32'). 3(c) shows the result of a 'smoothed' coarse-graining - using overlapping cells - 
which is more suitable for plotting graphs |3(d) is the same coarse-grained density as in Fig. |3(b)] from a different 



perspective; 3(e) is a close-up of a single coarse-graining cell, at the level of individual lattice points. The irregular 



nature of the underlying distribution at this level is clear; |3(f) is the same as Fig. 3(c) from a different perspective. 



For the 16-mode case, it was found in Ref. i22jthat H decays approximately exponentially as 



One of us (AV) has presented a theoretical estimate of the relaxation timescale r obtained by considering 



the behaviour of the second-time derivative of at t = (where H posesses a local maximum) [201 [2T] . 
As discussed in more detail below, it was shown that in the limit e — )■ 0, r scales inversely with e. Further 
estimates, or simply dimensional analysis, then suggested the rough formula [21] 

T ~ (9) 

Here e is the length of the coarse-graining cells and AE is the energy spread of the wave function. For 
reasonable values of e, this estimate was in rough agreement with the numerical value [22] . 

However, because H has a local maximum at t = 0, the estimate on the right-hand side of Eqn. [o]- obtained 
from the second time-derivative of if at t = - can only define a timescale that is valid close to t = 0. As 
we shall see, it cannot properly represent the timescale r associated with the subsequent (approximately) 
exponential decay. The results of this paper in fact show a scaling of r with Ai? that disagrees with Eqn. [9] 
- but which is in agreement with an improved estimate discussed below. 

The question of how r varies with the number of modes was not investigated by Valentini and Westman [22j 
owing to computational difficulties. That gap is filled in this paper. 



II. NUMERICAL SIMULATIONS 



In this work we compute the dependence of the relaxation time r on the coarse-graining length e and energy 
spread AE through explicit numerical simulations. An initially non-equilibrium probability density in a 2D 
infinite potential square well is evolved according to pilot-wave dynamics, using a wave function consisting 
of an out-of-phase superposition of the first M energy eigenstates (normal modes). For this choice of wave 
function, taking all the modes to have equal weight, we have AE ~ so we look for a dependence of the 
form 

r ~ e'^M'P (10) 
On the basis of Eqn. [9| for example, we would expect p = — 3 and q = —1. 

To study this system (and potentially others, since it is designed to be easily extendible) we have written 
a new computer code named 'LOUIS' [23] which uses pilot-wave dynamics to calculate particle trajectories. 
Given an initial probability density and wave function, LOUIS is able to use these trajectories to compute 
the probability density function and coarse-grained subquantum ii-function at later times. It is a ground-up 
reimplementation of the code used in Refs. [22] and [23] and is up to two orders of magnitude faster with 
significantly more capabilities. Currently it can treat infinite potential square wells in one, two, or three 



dimensions using a finite superposition of eigenstates to represent the wave function. The relative weights 
and phases of the eigenstates may be specified in the input or chosen randomly (but reproducibly, using 
preset seeds). The scale of the coarse-graining may also be set manually, outputting results for multiple 
coarse-graining lengths in a single run of the program. 

Since the results of interest involve calculation of the subquantum if- function, which in turn involves a 
numerical integration over the area of the potential well, the quantities in the integrand, p and |^| , must 
be evaluated on a regular lattice. In all calculations presented here, a 1024 x 1024 lattice is used covering a 
square two-dimensional cell of length tt. The 'coarse- graining length' e refers to the number of lattice points 
along one side of a coarse-graining cell. 

A. Details of the algorithm 

The LOUIS code uses de Broghe-Bohm trajectories to calculate how the particle probability density evolves 
from a given initial density. At each of a sequence of requested times, it evaluates the particle density and 
wave function at all points on the fine-grained lattice, and then applies coarse-graining on the requested scales. 
The coarse-grained iJ-function is calculated from these data at each timestep, and output files containing 
{t,H) pairs and {t,lnH) pairs are written to disk. The program calculates a straight-line fit using linear 
regression of the t vs. In H data; assuming exponential scaling, the gradient of this is the decay constant or 
relaxation time r. 

How do we calculate the density at a later time? We have seen that the ratio is preserved along 
trajectories, implying that the density at position x and time t may be calculated from 



where the positions Xq and x are points on the same trajectory, at times and t respectively. The value 
of 1^1^ can be calculated analytically at all positions and times, and p(xo,0) is a known function, therefore 
p(x, t) may be calculated directly once we know the trajectory endpoint x. This is the crucial relation used 
to calculate probability density functions from trajectories. 

In fact, certain practicalities require real calculations to be performed in a slightly different manner. The 
subquantum if-function is evaluated through numerical integration over the 2D box from a set of values of p 

and calculated at discrete points. Since accurate and efficient quadrature algorithms in few dimensions 
generally require the points to be sampled uniformly across the region, LOUIS starts with a uniform lattice 
at time t and exploits the time-reversibility of the dynamics to calculate particle trajectories backwards in 



p(x,t) 



*(x,t) 



2^(^0,0) 




|*(xo,0) 



time to t = 0. This ensures uniform samphng of p{^it) at t when the quadrature is to be performed. This 
has the unfortunate consequence that if p is required at a later time t' > t this 'backtracking' has to be done 
all the way to t = again: the data calculated at time t cannot be used again. 

The rate-limiting step of the LOUIS program is the numerical integration of the de Broglie guidance equation 
v(x, t) = VS'(x, t) (in atomic units) to compute the particle trajectories x(t). One may use a variety of 
standard algorithms; an excellent choice for these purposes is Runge-Kutta-Fehlberg [32] ■ Currently, the 
Schrodinger equation is not integrated numerically to compute the time-development of the wave function; 
instead, only finite superpositions of stationary states are used so the wave function can be evaluated exactly 
for any t. 

The velocity of the particle at any point may be computed from Im^^^ where the M-mode wave function 



As with all such algorithms, Runge-Kutta-Fehlberg basically involves adding small increments to a function - 
here x(t) - where the increments are given by derivatives (^ = v = VS* = Im^) multiplied by variable step 
sizes (here, a timestep At). In order to increase the accuracy, a tolerance is set for the maximum error on each 
step (the step tolerance); if the error is greater than this, a smaller timestep is used (subject to appropriate 
underflow checks). When the integration has been performed along the entire trajectory between the required 
initial and final times, the whole trajectory is recomputed with the step tolerance decreased by a factor of 
10. If the two final positions agree within a certain tolerance (the trajectory tolerance) , then the trajectory 
is kept. If not, the process is repeated with smaller and smaller step sizes until the trajectories converge, 
or until the step tolerance reaches a certain minimum value where the calculation will take too much time 
and the trajectory is flagged as failed. Failed trajectories are not used in the subsequent computation of the 
density. The proportion of failed trajectories rarely exceeds 1 in 1000 and their contribution to the overall 
error is negligible. In general, they are trajectories that come too close to a wave function node (where the 
velocity field diverges). 

Computational cost is the main limiting factor. The calculation of relaxation timescales is very computation- 
ally intensive, requiring many long, high-precision numerical integrations. Since the particle wave function 
and its gradient must be evaluated at each step in the integration of each trajectory, the complexity of the 



is given by 




(12) 




wave function contributes significantly to the time taken to perform the calculation; runs with larger num- 
bers of modes in the superposition are considerably more expensive. For example, a typical calculation on 
an elderly cluster of sixteen processors (seventeen evaluations of H between time and time 47r) took 9.6 
CPU-hours with a 9-mode wave function. A comparable calculation with 36 modes took 542 CPU-hours. 
The calculations we report here are for 4, 9, 16, 25, 36, 49, and 64 modes. Previous calculations were done 
exclusively with either 4 modes (Colin and Struyve [23]) or 16 modes (Valentini and Westman [22]). 



III. RESULTS AND DISCUSSION 



We begin by verifying that exponential decay is an appropriate model for the evolution of H under these 
conditions as was assumed in the definition of the relaxation timescale. Fig. |4] shows plots of In H vs. t 
for different coarse-graining lengths e, different numbers of modes M in the superposition, and different 
sets of initial phases Omn in the wave function (in order to obtain reproducible results the phases are fixed 
by a single 'preset' parameter in the LOUIS input file, which controls the seed for the random generation 
of a set of phases). In all cases we find a good straight-line fit to the data, validating the assumption of 
exponential decay in H over this range of conditions. This was not unexpected as exponential decay was 
previously demonstrated by Valentini and Westman [22] , though only for the case of M = 16 modes with a 
fixed coarse-graining length. 



A. Relaxation time r as a function of the number of modes M 

Fig. [5] shows plots of In r against In M for various coarse-graining lengths e where now logarithmic axes are 
used in order to search for a power law relationship of the form r oc M^. Error bars were calculated by 
running LOUIS six times for each (e, M) pair with different initial phases 9mn in each run; the mean of these 
was taken to be the best estimate for the timescale with the standard deviation taken to be the error bar. 

In the case of 4-mode simulations, the spread of values of r was so large that the error bar Ar could not 
be displayed on a logarithmic plot (hence the arrow at the base of the corresponding error bars in Fig. |5|. 
Some representative values of r were (980 ± 1600) for e = 4 and (100 ± 110) for e = 32. This large spread in 
the timescale for low M can be understood by considering the role of wave function nodes in the relaxation. 
The rapidly-varying velocity field in the vicinity of nodes is believed to be a significant driving force for this 
process |22j, and so the initial positions of the nodes are likely to affect the timescale (these initial positions 
are moved around when modifying the set of initial phases Omn)- The important point is that in larger 
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FIG. 4. Graphs showing that the model of exponential decay is robust under the conditions studied in this article. 
The horizontal axis on the first graph in |4(b)|is scaled differently to the others, to better fit the range of data. 



superpositions with larger M there are many nodes and the exact change in positions of nodes will have less 
effect because the average distribution will be similar. With a small superposition, maybe only containing 
one node, the initial position of this node will have a much larger effect on the subsequent relaxation, and 
so the six runs with different initial phases will tend to produce very different results. 

A full animation of the relaxation process can be seen online at Ref. [33]. The swirling vortices in the density 
surrounding the moving wave function nodes are quite striking, and one can obtain more of a visual sense of 
why the presence of nodes increases the chaotic nature of the trajectories. 

In Fig. |6j therefore, we choose to plot the same data as in Fig. |5]but with the points for M = 4 omitted, in 
an attempt to show the relationship between r and M in the regime where the distribution of nodes may 
be considered roughly fixed. Can we justify this? Arguably, points with such a large error will have little 
effect on the curve fitting, as points are weighted according to the inverse of their error. However, a better 
reason for neglecting the points for M = 4 is that the theoretical predictions are actually in terms of AE, the 
uncertainty in the energy of the wave function, rather than the number of modes, M, in the wave function. 
The approximation used in this analysis is that, /S.E ~ M^, which only applies for larger M. 

The best-fitting power law is shown on the graph in each case, and they are also summarised in Table |T) 
Clearly these results do not support the theoretical prediction (described as a 'crude estimate' in Ref. [22]) 
that the relaxation time should be proportional to . Instead they very strongly suggest a relationship 
of the form r oc M~^, to within the estimated error. This suggests that some of the approximations made in 
obtaining the prediction in Eqn. |9] are invalid. 

The relevant arguments used to obtain the apparently incorrect estimate of the scaling are set out in Ref. [21] . 
They begin by defining a relaxation timescale r in terms of the rate of decrease of H near t = via the 
following formula: 

1 1 f^m . ^^3^ 



r2 V dt2 



\ / 



The second derivative is used rather than the first derivative since =0 (^s may be shown analyti- 

cally pO]. The if-curve necessarily has a local maximum at t = 0. This property might, at first sight, seem 
to be incompatible with the observed exponential decay. But in fact, it must be the case that the exponential 
decay sets in soon after t = 0, and that in the limit t — )■ the decay is not exponential. The timescale of 



Eqn. 13 then applies only to this very short period immediately after t = 0. In other words, while Eqn. 13 
may well be a good estimate for the relaxation timescale in the limit t — t- 0, it cannot accurately estimate 
the time constant in the exponential tail - where the bulk of the relaxation takes place - and is therefore 
of little practical relevance. Another potential conflict in the derivation of Eq. |9] is the requirement for a 
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FIG. 5. Graphs showing Inr against InM, for varying coarse-graining lengths. Logarithmic axes have been used 
to identify a power law relationship, r oc M^. The error bars were estimated by running LOUIS six times for each 
(e, M) pair, using different initial phases. The values shown above correspond to the mean of these six runs, and 
the error bars are one standard deviation. The errors on the points for M = 4 cannot be properly represented on a 
logarithmic scale, as the lower bound is less than zero. There is reason to believe (see section [lTT| ) that these points 
may be excluded, and with such a large error their weight in fitting would be very small, so Fig. [6] shows the same 
results without these points. 
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FIG. 6. Graphs of r against M for all e, excluding the points for M = 4. The equations for the best power law fit 
(using the GNUplot implementation of the Marquardt-Levenberg nonlinear least-squares fitting algorithm) is shown 
on each of the graphs. |6(f) shows a comparison of the best fits for all e. It can be seen that for all e there is a 
roughly constant power law with index ~ — 1. 
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-1 06 ± 18 
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-1.09 ±0.17 
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-1.09 ±0.12 
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-1.08 ±0.03 


64 


-1.05 ±0.03 



TABLE I. Summary of results for a relationship between r and 
M, for various coarse-graining lengths. The error is estimated 
from the straight line fit, which takes into account the errors on 
the points used for the fit. The values are clearly not compatible 
with r oc as in Eqn. [oj but are compatible with r oc M^^, 

within error, supporting the relation in Eqn. 14 





(a)fine- grained 



(b)coarse-grained 



FIG. 7. Comparison of the fine-grained 1^*1 and the coarse-grained approximation, demonstrating that there may 
be significant variation over the length of coarse-graining cells. 



coarse-graining length e so small that the velocity field varies little over the length of a coarse-graining cell. 
The derivation works by considering the dependence of r on e in the limit where e — )■ 0, then applying 



dimensional analysis to find the other dependencies. The inverse dependence of r (as defined by Eqn. 13) 
on 6 must hold in this limit, since we have an analytic proof of it. However, unfortunately, the limit seems 
too restrictive to be of practical use. In the cases studied here the wave function apparently varies rapidly 
enough on that scale - particularly with larger numbers of modes in the superposition - to ensure that we 
are not working in this limit at all. For example. Fig. [T] demonstrates the effect of coarse-graining {e = 32) 
on one of our 64-mode wave functions. This shows the magnitude of rather than the velocity field, but 
the length scale over which the latter varies should be even smaller than the length scale over which the 
former varies. We may thus conclude that the velocity field is likely to vary significantly over the length 
of one coarse-graining cell, at least under some of the conditions studied here. Taking all these facts into 
consideration, it is not surprising that the results do not confirm Eqn. |9| whose domain of validity is probably 
simply too narrow to be of practical use. 



We now provide a theoretical justification for tlie observed relation r oc M~^. Let qr{r = 1,2) denote x or 
y. We shall proceed by considering an upper bound on the (equilibrium) mean displacement 6qr of particles 
over an arbitrary time interval [ti,tf]. A relaxation time may then be defined, in the case of the infinite 
potential square well, by the condition that relaxation will occur over timescales r such that the said upper 
bound becomes of the order of the width L of the potential well. As we shall now show, the timescale will 
then be 



where E is the mean energy, and E oc M^. 

There now follows a derivation of the upper bound on the mean displacement 6qr. The derivation is based 
on Ref. |3H but with some differences that are highlighted below. 

First, note that the final displacement 6qr (t/) has modulus 




(14) 




(15) 



where (jr is the component of the de Broglie-Bohm velocity in the direction. 



The equilibrium mean {\6qr (i/)|)eg then satisfies 




(16) 



The equilibrium mean speed, {\qr (^)|)e„ is 




(17) 



Using the fact that, for any x, (x) < >y (x^), we have 




(18) 



Now note that, starting from the guidance equation: 




(19) 



where tt^ is the momentum operator conjugate to qr, and {nr) denotes the usual quantum expectation value 
for the operator vr^. The last equality follows from the relations 
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dx dy^— — 
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(20) 



and 



dqr dqr \ dqr J \dqr J 



W U- ■ (21) 



Since (9 |^| /dqr)^ > 0, we then have 

^'(l^^Oe (22) 



m 



1 r'f 



and so 

{\^1rm)e,<-l dt\l i^'r}- (23) 

We also have 

{^tt)^ < 2m (24) 
where Wr denotes the x- or y- part of the mean Hamiltonian, with Wr cx M^. Hence, 

(I'^^r (t/)|)e, < [ ' dt-V^M^. (25) 

Since W,. is time-independent, and setting = and tf = t we have 



(l'^^^(^)l)e,<n/^- (26) 



Setting the right-hand side to be of order L, and noting that E ^ Wr, then indeed yields the relaxation time 
in Eqn. [Ml whose inverse scaling with M is in agreement with the numerical results presented above. 



As mentioned, this derivation is based on that in Ref. |34| but differs in some respects. The purpose of the 
analysis in Ref. |3l] was to derive a condition for the suppression of relaxation in expanding space (here we 
are only concerned with static space) and the condition for relaxation was that the mean displacement Sqr - 
for field degrees of freedom in Fourier space - should be comparable to (or greater than) the quantum spread 
Aq^. In the analysis above, the only degree of freedom considered is the spatial displacement of a particle 
in the potential well, the constraints of which slightly change the condition for relaxation. Regardless of the 
spread in the wave function the particle cannot move beyond the confines of the well, so the condition used 
for relaxation is that the mean displacement of a particle is comparable to (or greater than) the size of the 
potential well. 
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TABLE II. Summary of results for a relationship between r and 
e for various numbers of modes, shown graphically in Fig. [8j 
The errors are estimated in the same way as the results in Ta- 
ble |lj These values are not compatible with Eqn. [oj r oc e^^, 
nor do they appear to be compatible with a consistent power 
law. 



B. Relaxation time as a function of coarse-graining length 

In Fig. [8} the relaxation time r is plotted against coarse-graining length e for various M, again with log- 
arithmic axes and with the data for M = 4 excluded for the same reason as before. The best-fit power 
law is shown on the graphs, and the results are summarized in Table [Tij It is evident from these data that, 
although the relaxation time evidently decreases with increasing e, the data do not support Eqn. [9] nor are 
they particularly suggestive of a constant power law. A weak dependence of order r ~ e^^^'^ is observed. 



Eqn. M predicts a power law relation t cc e ^ while Eqn. 14 predicts no dependence. 



As was the case with the M-dependence of r, the dependence on e is not in fact expected to take the form 
in Eqn. |9} for two reasons. First, the derivation of Eqn. |9]is based on a definition of 'relaxation time' that 
applies only when t — )■ 0; it does not apply to the exponential tail of H, where most of the relaxation takes 
place. Second, the velocity field can vary significantly over a coarse-graining cell, contrary to the assumption 
made in the derivation of Eqn. [9j Indeed the graphs in Fig. [8] appear to show a systematic concave character 
rather than a straight line, apparently suggesting a small systematic deviation from a power law model. The 
concave curvature is at least consistent with a power law r oc e'^ in the limit e (which must hold 
as we have noted), since at some point to the left of the data in any of the figures the gradient ought to 
approach (or pass through) —1. A systematic study at smaller coarse-graining lengths would however be 
rather difficult in terms of CPU time because of the need for a significantly finer lattice and much greater 
overall number of lattice points. 



The lack of a dependence on e in Eqn. 14 is not surprising since the analysis in section III A and [31] uses 
a different definition of the timescale. This definition does not consider the ff-function nor is there any 
necessity to even mention coarse-graining. The weak dependence observed in the numerical simulations 
should be interpreted as an effect outside the scope of this prediction rather than one in conflict with it. 
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FIG. 8. Graphs showing Inr against Ine. As in Figs. [5]and[6j a logarithmic axis has been used, so as to identify a 
power law relationship. Also, as in Fig. [6| the datapoints with M = 4 are excluded. The errors are estimated in the 
same way as those in Fig. [g} It is clear that the power law does not conform to Eqn. [oj r oc e^^, nor does there seem 
to be a consistent power law. The straight line fit is not as convincing in this case as in the t vs M graphs, except 
for the case of M = 9. 



IV. CONCLUSIONS 



The numerical simulations performed in this work demonstrate clearly and unequivocally the tendency for 
Born-rule distributions to arise spontaneously as a consequence of ordinary pilot-wave dynamics, even for a 
system as simple as the electron in a two-dimensional potential well. Contrary to popular belief therefore, 
the Born rule does not have to be introduced as a postulate of non-relativistic quantum mechanics. What 
is the price paid for this? We must suppose merely that particles have well-defined positions (and hence 
trajectories) continuously rather than only when a position measurement is performed. 

The main technical result of this work is the emergence of the relationship r oc showing the dependence 
of the relaxation time on the number of modes in a superposition. This result seems fairly robust under 
the conditions studied here. In general terms, the faster relaxation times for larger M are due to the 
greater number of free-moving nodes in the pilot wave which act as a source of vorticity and increase the 
chaotic nature of trajectories. Our numerical result for the scaling conflicts with the previous theoretical 
prediction, Eqn. [9} but agrees with an alternative theoretical analysis presented here in section III A As we 



have discussed, the assumptions made in deriving Eqn. |9]were probably too restrictive for it to be of practical 
use. In particular, the defined timescale is relevant only close to t = 0, and does not apply to the exponential 
tail of the H function, where most of the relaxation takes place. 

Our simulations reveal no well-defined scaling for the relaxation time as a function of coarse-graining length 
r {e), other than a possible weak-dependence of the order of r ~ e~^^^. This also differs from Eqn. [o] and this 
is probably simply because the coarse-graining cells are too large for the derivation of Eqn. |9] to be valid. It 
is possible that by decreasing e a behaviour conforming better to Eqn. [9]( r oc e^^) could be observed, and 
it would be interesting to see at what length scale this begins to emerge. 

Physically speaking our results suggest very short relaxation times with a range of values observed for 
r between about 1 and 1000. Using natural units c = h = 1 and an electron with mass m = rrie = 
1 this corresponds to relaxation times of the order of 10^^^-10^^'^s. This is consistent with our current 
understanding of quantum mechanics and modern experimental investigations in which no deviation from 
quantum equilibrium is observed. If the initial state of the universe corresponded to a non-equilibrium state, 
one might then assume that almost all deviations from the Born rule will have been quickly washed out. 
However, it may be (see Refs. HHI and [M] and works cited therein) that relic non-equilibrium from the early 
universe could be observed today, either directly or by its imprint on the cosmic microwave background. In 
the future we intend to modify the LOUIS program to deal with more realistic wave functions, multi-particle 
systems, and expanding spaces with the intention of improving some of the predictions outlined in Refs. 



and 1311 More precise predictions for the effect on the CMB of non-equihbrium in the early universe could 
lead to experimental tests of the de Broglie-Bohm formulation of quantum mechanics. 

How might these results generalize to more complex systems? In the original version of the subquantum 
iZ-theorem published by one of us (AV) in 1991 [16], relaxation was considered for a theoretical ensemble of 
complex A'^-body systems. Once equilibrium is reached for such an ensemble, it can be (and was) shown that 
extracting a single particle from each system resulted in single-particle sub-ensembles that obey the Born 
rule. The original view expressed in Ref. dni was that, realistically, relaxation would take place efficiently 
only for many-body systems, and that the Born rule for single particles would be derived by considering how 
these are extracted from more complex systems. But in fact, in practice, there is an efficient relaxation even 
in simple two-dimensional one-electron systems, and there appears to be no need to appeal to a complex 
iV-body 'parent system'. 

Note that in a strict account of our world, say in the early universe, it could well be that all degrees of 
freedom are entangled, so that there is no actual ensemble of independent subsystems with the same wave 
function. However, one can still talk about a theoretical ensemble of universes, each with the same universal 
wave function, and consider the evolution of its distribution. (One could also consider a mixed ensemble of 
universes, and apply our discussion to each pure sub-ensemble.) 

It is sometimes suggested that it is problematic to consider probabilities for the 'whole universe'. And 
yet, cosmologists are currently testing primordial probabilities experimentally by measuring temperature 
anisotropics in the cosmic microwave background. By making statistical assumptions about a theoretical 
'ensemble of universes', cosmologists are able to test probabilities in the early universe, such as those predicted 
by quantum field theory for vacuum fluctuations during inflation. (For a detailed discussion of this in a de 
Broglie-Bohm context, see Ref. |19j|.) One can question what the ensemble of universes refers to. Is it a 
subjective probability distribution? Or, is the universe we see in fact a member of a huge and perhaps 
inflnite ensemble, as is the case in theories of eternal inflation? Those are interesting questions, but only 
tangentially related to the ongoing experimental tests. It is also important to bear in mind that there is 
much that is not known about cosmology, so the treatment should be kept independent of cosmological 
details as far as possible. What we do know is that all the particles we see today are, or were, in complex 
superpositions (whether entangled with other particles or not), and it appears clear from the simulations 
that such superposition yields rapid relaxation - if it is rapid in two dimensions, one would expect it to be 



even more rapid for 3A^ dimensions. 
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